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SIMILARITIES AMD DIFFERENCES IN SIGNAL PROCESSING 
FOR RADAR AND SONAR 

David Creasey* 



1 INTRODUCTION 

1.1 Information 

Both sonar and radar systems are types of communication 
systems. Shannon states that in a single ideal communication 
channel, the maximum information output is given by [1] 

I < 2BTlog 2 ( 1 + s/n ) bit (1) 

where B = the bandwidth of the channel, 

T = the observation time and 
s/n = the signal-to-noise ratio. 

Neither radar nor sonar are ideal information channels, but 
Equation 1 does indicate what we must do in order to get more 
information from the system. 

Firstly, note that if the signal-to-noise ratio approaches 
zero then the information that can be derived from the system is 
also virtually zero, (log 2 l - 0). Put another way, if you put 
rubbish in you get rubbish out. The system designer and the 
system user must always ensure that there is signal to process. 

Secondly, a long observation time will produce more 
information than a short observation time. Thus if we can 
integrate the output of a system over along time then we shall 
obtain some processing gain. 

Thirdly, always use a signal with a wide bandwidth. In some 
ways, this contradicts the requirement for a high signal-to-noise 
ratio since noise increases as the bandwidth increases. However, 
the signal-to-noise term is included in a logarithm, so there is 
a gain in an information sense by using a high bandwidth. 

Often radar and sonar systems use an array of elements. If 
the array elements are sufficiently spaced, more than half a 
wavelength apart, then each element can be considered to form a 
single channel. This means in a system with an array of M 
elements, each suitably spaced, there are M separate channels. 
Shannon's simple information limit can then be increased by a 
factor of M to give 

I < 2BTMlog 2 ( 1 + s/n ) bit (2) 

1.2 Active and Passive Systems 
1.2.1 Active systems. 

If energy is propagated into a medium, an object within that 
medium will intercept and re-radiate some of the transmitted 
energy. A receiver will then be able to detect the presence of 
the object by observing the reflected energy. Such systems are 
known as echo-ranging systems since they measure the two-way 
propagation time and, from a knowledge of the speed of 
propagation, the range of the target can be calculated. 

The most obvious difference between a radar and a sonar 
system is the difference between the speed of electromagetic 
propagation for a radar (300 x 10 6 m/s) compared to the speed of 
sound in water (1.5 x 10 3 m/s). For a target at a range of 
15 km the two-way propagation time for a radar is only 100 ^is. 
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For a sonar target at the same range the two-way 
propagation time is 20 s. Things happen very much more 
slowly in sonar than they do in a radar. 

There are three sources of interference in both radar and 
sonar systems. Firstly, there is the background noise. A signal 
that is always present due to many causes. These causes can be 
naturally occuring noise, or, in the other extreme, it may be 
noise deliberately generated to confuse the echo-ranging system. 
The system must ensure that the output of the receiver can have a 
sufficiently large signal component so that it may be detected in 
the presence of noise. This can often be accomplished by 
increasing the quantity of energy radiated by the transmitter. 

There are often many scatterers both within the medium and 
at its boundaries. These scatterers are often individually small 
but there are many of them so that they can combine to produce 
large but very variable signals at the receiver. Such signals are 
called clutter in radar and reverberation in sonar. Increasing 
the transmitted signal level does not improve the 
s igna 1-to-r ever berat ion level since they are both proportional to 
the level of the signal transmitted. A large amount of 
directionality in the receiver or the differences in Doppler 
shift between the signals backscattered from the target and the 
clutter /reverberation can be used to differentiate between 
signal and clutter /reverberation. 

The third form of interference signal is very similar to 
clutter /reverberation. Energy can often propagate over different 
paths, for example by a direct path and by reflection from a 
boundary. Differences in propagation time can be used to 
differentiate between the direct path and the reflected path. 

The signal-to-noise ratio at the receiver can be evaluated 
from a number of factors. These result in an energy balance 
equation where the level of the echo signal is compared to the 
interference (noise) . The resulting radar equation and sonar 
equation are both derived by taking the transmitted signal level 
multiplying by the losses expected in the two-way transmission 
path and further multiplying by the ratio of the energy reflected 
from the target to that intercepted by the target. This gives the 
level of the echo signal at the receiver. This is then divided by 
the level of the noise to give the signal-to-noise ratio. Often 
the series of multiplications and divisions provides a very 
cumbersome representation. Alternatively, the ratio can be 
expressd in a logarithmic form and usually this is done by taking 
the terms individually and expressing them in logarithmic form. 
This results in an equation of the type 

S/N = SL - 2 PL + TS - NL dB (3) 

Equation 3 appears very simple but it hides many components. For 
example the source level, SL, will depend upon the directivity of 
the transmitter; the propagation loss, PL, will contain 
components that comprise spreading losses and absorption losses; 
the target strength, TS, is aspect dependent for many practical 
targets; and the noise level, NL, will depend upon the noise 
spectral density, the bandwidth and possibly the directional 
properties of the receiver. In radar and sonar the basic 
formulation of Equation 3 is identical but the detailed 
components within each of the terms differ . 

Active sonars operate over frequencies from a few hundred 
hertz to about 1 MHz. Radar on the other hand operates from r.f. 
frequencies of say 10 MHz up to frequencies just below the 
infra-red region at about 30 GHz. Thus a second major difference 
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between radar and sonar is the operating frequency. Strangely the 
wavelengths used overlap because of the major difference of 
propagation speed. The relationship between velocity, c, 
frequency, f, and wavelength, X, is 

X = c/f (4) 

Thus, the wavelengths for an active sonar range from about 3 m 
down to 1.5 mm. The range of radar wavelengths are from 30 m 
down to 10 mm or so. 

1.2.2 Passive operations 

The need for covert operations underwater has lead to the 
development of sophisticated systems that rely upon the noise 
created at the target by its propellers, on-board rotating and 
reciprocating machinery and flow noise. 

Radio telescopes are effectively passive radar systems 
although they are not often viewed as such. The nearest direct 
passive equivalent usually associated with radar area electronic 
support measures (ESM) where listening equipment detects 
radiation from active radars and communications systems. Sonar 
intercept equipments are of course used in the underwater 
environment . 

The operating frequencies of ESM and sonar intercept 
equipments are dictated by the frequencies used in active radars 
and active sonars. Ideally, passive sonars would operate at 
frequencies from 10 Hz or so and up to about 3 kHz so as to 
encompass all the useful spectrum of the radiated signatures of 
passive targets. The difficulty of operating at the frequencies 
of around 10 Hz is the high ambient sea noise and the difficulty 
of operating with large apertures so as to obtain directional 
information. 

The major interference component in a passive sonar is the 
ambient noise of the sea. Multi-path propagation can cause 
problems when narrow-band operation results from such operations 
as spectrum analysis. Here the signals from different paths can 
combine either constructively of destructively and the signal 
levels perceived by the sonar can fluctuate wildly. However, the 
basic premise of a passive sonar is often wideband operation so 
that over the wideband the signal level is more constant. 

1.3 Objectives 

The objectives of radar and sonar are threefold. They must 
detect, locate and classify targets so that effective 
actions, such as the derivation of a fire-control solution, can be 
obtained. All of these objectives require high signal-to-noise 
ratios to be effective. 

Detection requires that the target signal should be well 
above any interference. If that interference is clutter or 
reverberation, it may be necessary to look for differences in 
Doppler between the target and the interference. When an active 
system is noise limited and it may be impossible to increase the 
peak transmitted signal level. One solution is to code the 
transmitted signal, for example, changing the frequency during 
the pulse. On reception the echo signal plus noise is cross 
correlated agianst a replica of the transmitted signal. The echo 
should produce a high degree of correlation while the noise 
should be poorly correlated. 

The range, bearing and heading are the three basic 
parameters that an echo ranging system can produce. The range is 
obtained simply from the time between pulse transmission and echo 
reception. Bearing is obtained from the directional response of 
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the system arrays. Heading information is obtained from the 
changes in range and bearing observed over a period of time. 
Doppler shift will give radial speed directly. 

In a passive system, range is not directly available and must be 
computed from a number of sets of data samples either separated 
in time or space. 

Classification is a very sensitive topic. Suffice to say, if 
a very short pulse is used the amplitude highlights in the echo 
signal may give clues about the target shape. Alternatively the 
behaviour of the target observed over a period of time may be the 
only clue available in an active system. In passive sonar, the 
received signal spectrum will contain a number of discrete lines 
that may be used to classify the target. 

Thus radars and sonars both need to generate signals in the 
receiver above interference signals. They need to look at the 
spectrum of the received signals. Time integration is used in 
replica correlators and both need to process the data received 
from arrays of elements. 

1.4 Bandwidth 

The speed of propagation causes the major difference between 
the digital processing of radar and sonar signals. This speed 
difference is directly responsible for the differences in system 
bandwidth. The bandwidth in an active system is fixed by the need 
to obtain a good range resolution. The range-resolution cell is 
equal 0.5cT where T is the effective pulse length and c is the 
velocity of propagation. The system is matched when it has a 
bandwidth, B, just sufficiently wide to accommodate the pulse 
length. This is usually expressed in the 
simple equation 

B = 1/T (5) 

Table 1.1 gives comparative parameters for typical radars and 
active sonars. 



Tft P X iE i ,i 



Activitv 




Operating 

frecruencv 


Pulse 

length 


Bandwidth 


Range res'n 
cell 


Long Range: 


Radar 


1.25 GHz 


1 US 


1 MHz 


150 m 


Sonar 


4 kHz 


1 ms 


1 kHz 


0.75 m 


Weapons System 


Radar 


15 GHZ 


30 ns 


33 MHz 


4.5 m 


Sonar 


40 kHz 


100 us 


10 kHz 


75 mm 



This bring us to the question about how signals should be 
sampled prior to quantisation and being processed digitally. It 
is possible for signals from an active soanr and a passive sonar 
to be sampled directly as baseband signals at a rate something 
slightly above the Nyquist rate. For the long range sonar 
operating over the band 3.5 kHz to 4.5 kHz then a sampling 
frequency of about 10 kHz will suffice. Even the signals from a 
weapons sonar could be sampled as a baseband signal at about 
sampling rate of 100 kHz. 

However, analog-to-digital converters operating in the 
gigahertz region are very rare. To accommodate radar signals the 
signals have to be modulated down to baseband frequencies. When 
this is done, it is usual to modulate down to baseband using a 
local oscillator working at the centre frequency of the radar 
signal. The signal then needs to be modulated with both the sine 
and the cosine of the carrier in separate channels to produce 
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in-phase (I) and quadrature (Q) signals. The sampling frequency 
in each channel can then be halved but as there are two channels 
the number of samples is still equal to those rquired by the 
Nyquist sampling theorem. Table 1.2 compares sampling rates for 
radar and sonar signals. 



TABLE 1.2 



Activitv 




Operating 

freauency 


Bandwidth 


Min. Sampling 
Freauencv 


Long Range 


Radar 


1.25 GHz 


1 MHz 


1 MHz in both 
I & 0 channels 


Sonar 


4 kHz 


1 kHz 


2 kHz as base- 
band signal or 
1 kHz in both 
I & 0 channels 


Weapons System 


Radar 


15 GHZ 


33 MHz 


33 MHz in both 
I & 0 channels 


Sonar 


40 kHz 


10 kHz 


20 kHz as base- 
band signal or 
10 kHz in both 
I & 0 channels 



Note that for many situations sonar signals need to be sampled in 
quadrature and apart from the sampling frequencies involved radar 
and sonar signals are often treated similarly. 

2 SPECTRUM ANALYSIS 

2.1 System Model 

Spectrum analysis is one of the basic processes used in 
signal analysis. The technique is based upon a linear-system 
model that assumes all signals consist of a summation of cosine 
waveforms of different frequencies. So that time delays can be 
incorporated each cosine term has an associated phase shift. 
Mathematically this can be written simply as 

x(t) = IX n cos(w n t + <ki) (6) 

Alternatively, the cosine term can be expanded to give 

x (t) = KX n cos4>n.cos(0 n t - X n sin<J>n. sinw n t) ....(7) 

= £a n cosco n t - Ibnsinwnt .... — (8) 

Complex waveforms, such as noise, can consist of an infinite 
number of sines and cosines. The problem in spectrum analysis is 
to determine the terms X n and 4>n or the alternative pair a n and 
bn- We will now discuss how a practical spectrum analyser could 
be constructed. 

2.2 Spectrum Analysers 

The simplest method of realisation would be to feed the 
signals through a set of band pass filters each with a narrow 
pass band and centred at different frequencies, see Figure 2.1. 
This has three disadvantages: 

i) cost 

ii) the required phase term is not available 

iii) such filters may be impossible to design. 

The system could be made simpler by mixing the signals with 
a local oscillator and selecting one of the sidebands in a single 
filter, see Figure 2.2. This has the advantage that only one 
filter is required. To check all frequencies the local oscillator 
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Fig. 2.1 




Fig. 2.2 

frequency would need to be altered. This is how a radio or 
television receiver works. 

The design of the band pass filter could still be a problem 
if the centre frequency is too high. A low pass filter could be 
substituted for the band pass filter if the local oscillator 
frequency is made equal to the frequency of the signal being 
measured, see Figure 2.3. 




Fig. 2.3 

If the local oscillator is say 2.cosio n t, and the input 
signal is a n cosio n t - bfcSinWnt, the modulation process produces a 
waveform 

a n + a n cos2to n nt - bnsin2w n t 

The low pass filter removes the last two terms (the upper 
sidebands) and only passes a n (the lower sideband) . The problem 
with this realisation is that there is no measure of fc^. This may 
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be overcome by incorporating a second channel where the local 
oscillator is 2.sinco n t, see Figure 2.4. There are now two outputs 
one from each channel, a n and b n . 




2.3 Mathematical Representations 

To represent this process mathematically, we need to express 
the modulation processes as simple multipliers and the low pass 
filters as integrators. The two channels also need separating by 
using a multiplier, say a, in front of the sine channel. Hence, 
for the time waveform x(t), the component, X(co) , with angular 
frequency co can be written 

X(io) = Jx(t) . (cosiot + a.sincot).dt (9) 

-•© 

Note that the factor 2 in the multiplier has not been used. This 
is because in Eqn 7, there are two outputs, one at an angular 
frequency -co as well as the expected component at +co. 

2.4 The Fourier Integral Equation 

The a indicates that two channels are being used and that 
the resulting two terms are not added directly. The two channels 
act in quadrature to one and other. Of course, engineers and 
mathematicians commonly use "j" or "i", the square root of -l to 



represent quadrature components . In fact we can write a = - j so 
that Eqn. 9 becomes 

X (co) = Jx(t) . (cosiot - j.siniot).dt (10) 

-A© 

Remembering that by de Moivre's theorem 



cosiot - j.s incot = e^ 0 ^ 
a more compact form of Eqn. 10 becomes 

AO 

X(u) = Jx(t) .e"^ wt .dt (11) 

— AO 

This is the Fourier integral equation, one of a pair of equations 
that form the Fourier transform pair [2]. The second equation 
forming the pair is 
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X(t) = u/2lt).f X(w) .e^.dw (12) 

Eqn. 11 enables a time waveform to be represented by a collection 
of sine and cosine waveforms and produces the so-called spectrum. 
Eqn.12, on the other hand provides the reverse operation where 
the spectrum can be converted back into the time waveform. 

2.5 The Discrete Fourier Transform (DFT) Pair 

To calculate the Fourier transform in a digital form 
requires that the input signal is a sampled set of data. Thus 
x(t) becomes x(k.At) , where samples are taken at instants of time 
At apart and for brevity we write 

x(t) = x(kAt) = x k 



2.5.1 Band limiting 

If we assume that we are dealing with a baseband signal with 
bandwidth B and sampled at the Nyquist rate, the time samples are 
taken at instants At = 1/2B. Thus in Eqn. 11 and Eqn. 12 we can 
replace the time term t in the exponential by k/2B. 

2.5.2 Time limiting 

Any set of data being processed by a computer must be 
limited in length. Time signals found in radar and sonar have to 
be time limited. Suppose that we place a time window around our 
time waveform so that 

x (t ) = 0 for t < -0.5T, 

x(t) = Xk for -0.5T < t < +0.5T, and 

x ( t ) = 0 for t > 0.5T. 

The Fourier transform is a linear process so that superposition 
applies. This means that we can apply one component of signal at 
a time to Eqn. 11. The result using this signal can be evaluated 
and then added to similar results obtained by applying other 
signal components to Eqn. 11. Hence, suppose that we consider 
that within the time window -0.5T to +0.5T 

x(t) = Aicoswjt 

Applying this to Eqn. 11 and evaluating with the integral limits 
+/-0.5T gives the standard result 

X(<o) = 0.5A!T. [sintO.Sfco-WiJTj/fO.Stw-WiJT) + 

sinfO.Sfw+WiJTj/fO.Bfw+CDiJT}] (13) 

Eqn. 13 is sketched in Fig. 2.5. It can be seen that the 
function X(io) has two peaks, one at w = -«! and the second at co 
= +<*)!• This gives rise to the concept of negative frequency which 
we shall deal with in a moment. However, if we consider Eqn. 13 
more closely we note that the positve frequency component peaks 
at io= +( 0 i and that this component first passes through zero 
when the sine function in the sinX/X term equals +/-*. Thus zeros 
occur at at to = +t 0 i - (2 tt/T) and at to = +to x + (2ir/T) . The size of 
the frequency resolution cell can be defined as 1/T Hz. 

Remembering the 2n multiplier involved between frequency and 
angular frequency, the resolution cell size is equivalent to 
placing sinX/X functions along the frequency axis so that the 
peak of one function coincides with zeros for the other 
functions, see Fig. 2. 6. 

Returning now to Eqn. 11 and Eqn. 12, the continuous 
spectrum X(to) must be replaced by samples of that spectrum at 
frequency spacing 27r/T radian. These samples we shall call 
Fourier coefficients, A r . In Eqn. 11 and Eqn. 12 the to in the 
exponential terms must also be represented by samples. Hence, we 
replace w by 2irr/T. 
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0.5Aj T 



0 

ANGULAR FREQUENCY 



Fig. 2.5 



l/T 





x(cj) 




2.5.3 The number of samples. 

It is impossible [2] to limit the frequency bandwidth of a 
signal while simultaneously time limiting that same signal. A 
time-limited signal requires an infinite bandwidth and vice 
versa. However, if the number of samples taken is sufficiently 
large the resulting errors are sufficiently small for the 
equations to provide a good approximation. Combining the ideas of 
Section 2.5.1 and Section 2.5.2, we see that the time waveform of 
duration T has been sampled at instants separated by 1/2B. Thus 
the number of time samples is 2BT. If the one-sided bandwidth is 
B we also note that there is a band of positive signals and a 
band of negative signals. This means that the total two-sided 
bandwidth is 2B. Samples of the spectrum are taken l/T apart so 
the total number of Fourier coefficients is also 2BT. This 
indicates that the total information in the system remains 
unaltered by the transformation. 

The integral equations that define the Fourier transform 
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pairs must be replaced by summations in the discrete format. If 
we write 2BT = N, the total number of data samples in the time 
and frequency domains, we obtain equations of the form 



A r = 1/N N E 1 x k .e -2 * jrk/N (14) 

k=0 

x k = N £" Al .e +2 *^/N (is) 

r=0 



These two equations represent the discrete Fourier pair [2]. 

The 1/N multiplier can be found in the equation for x^ 
rather than in front of the equation for A r . Alternatively, some 
authors use (1/N) 0 * 5 [3]. The reason for using Eqn. 14 and Eqn. 

15 is that Eqn. 14 does at least have a physical significance. 

For example, with r = 0, Aq is simply the average value. Thus the 
zero- frequency term is equivalent to the d.c. value. For other 
values of r, the time samples, x^, are multiplied by samples of 
cosine and sine in orthogonal channels. If a frequency component 
equivalent to r is present, the summation and divide by N, then 
averages the lower sideband or base-banded component. Hopefully, 
the upper sideband components resulting from the multiplication 
will average to zero. 

2.5.4 Relationship to the Fourier series. 

If a continuous signal, x(t) , is periodic with a period 
27 t/w 0 , the signal may be expressed as a complex Fourier series [2] 

jrco 0 t 

x (t ) = Ic r ,e (16) 

r=-o o 

Comparing Eqn. 15 and Eqn. 16, it can be seen that the two 
equations are almost identical. In Eqn. 15, the time waveform is 
continuous whereas in Eqn. 14 the time wavefrom is sampled and 
time is represented by the variable k/2B. The fundamental 
frequency in Eqn. 16 is <o 0 /2 it. The equivalent fundamental 
frequency in Eqn. 15 is r/T. 

It is important to note that the discrete Fourier transform 
is simply a band-limited discrete version of the complex Fourier 
series. Such a representation requires that the signal being 
analysed is a PERIODIC signal with a period T. 

In radar and sonar it is very unlikely that the signals 
being analysed will be periodic and as such the DFT is only an 
approximation. The non-periodic nature of the signals analysed 
gives rise to errors in the estimation of the Fourier 
coefficients. In particular, the Fourier coeficients are spaced 
along the frequency axis at spacings of 1/T Hz. If a component 
does not fall at one of these discrete frequencies, r/T, it will 
cause all the frequency bins (all values of r) to have an output. 
These outputs will be modulated in amplitude by the sinX/X 
envelope. Thus if the frequency component is at (r+5) /T Hz, where 
5 is a positve fraction less than unity, bins r and r+1 will have 
major outputs from the forward DFT operation. Bins r-1 and r+2 
will also have significant outputs and these might be only about 
-13 dB below the expected output. 

2 . 6 Negative Frequencies 

The concept of negative frequency is sometimes difficult to 
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comprehend. It may be explained pictorially by reference to 
Fig. 2.7. In Fig. 2.7(a), the arrow represents a vector of unit 
amplitude with a phase 4>. If a perpendicular is dropped to the x 
axis the distance along the x axis is cos<|). A second line drawn 
to the y axis will be offset along the y axis by an amount sin<t>. 
If 4>-cot the vector becomes a phasor that rotates at an angular 
frequency co. In representing cosiot as a phasor, the problem is 
to remove the sintot component on the quadrature axis. This is 
done by having two phasors each of amplitude 0.5 rotating in 
opposite directions, see Fig. 2.7(b). Along the x axis the 
phasors add vectorially to give coswt. Along the y axis the 
resulting components add to give zero. 





+« t 



-u t 



(a) 



(b) 



Fig. 2.7 

Effectively, Fig. 2.7(b) represents the function 

0.5[e3 wt + e-3 wt ]. 

This is of course the expression for cos$ derived from de 
Moivre's theorem. 

A similar phasor diagram will produce a representation for 
sin<|>. From this second phasor diagram sin$ can be seen to be 
equal to 

[e jMt - e" jwt ]/2j (17). 

2.7 Doppler Effects 

In active systems, relative movement between the the source 
and receiver on one hand and the target on the other will produce 
a Doppler shift. The received pulse will be lengthened or 
shortened relative to the transmitted pulse. This in turn 
produces a change in the received signal spectrum compared to the 
spectrum of the transmitted pulse. For a relative radial velocity 
u (considered to be positive for a closing target) , the 
fractional change in angular frequency is given by 

A(d/co 0 = 2u/c (18) 

where c is the velocity of energy propagation. 

Table 2 . 1 shows the Doppler shifts expected for a relative 
radial velocity of Im/s in typical sonar and radar systems. 

Radial velocities in sonar can be from zero to 40 m/s. In 
radars the radial velocities can be from zero to values in excess 
of 600 m/s. Thus, the Doppler shifts in radar and sonar are 
similar in magnitude. 
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TABLE 2.1 



System 


Carrier frequency 
(Hz) 


Doppler shift 
per lm/s radial 
velocity (Hz) 


Sonar 


3 X 103 


r 4 




30 X 103 


40 


Radar 


j~ 1.5 XlO 9 


0.1 




15 X 10 9 


1 



It is important to note that differences arise when the 
spectrum of the signals are considered. Eqn. 13 can be adapted 
for this purpose. Consider a pulse of a carrier signal of angular 
frequency coq and of pulse width T. Eqn. 13 is then the spectrum 
of this pulse, see Fig. 2.5. However, in sonar and radar systems, 
the pulse itself is repeated at a regular pulse repetition 
frequency. This pulse repetition frequency has to be small 
enoegh to allow the energy to propagate to maximum range and be 
reflected from targets at that range. Thus the pulse repetion 
period, T p , for a system operating to a maximum range R is given 
by 

T p = 2R/c ( 19 ) 

The periodic pulsing of the transmit signal produces a line 
spectrum with a spacing between the lines equal to 1/T p . The 
spectrum of Fig. 2.8 shows the positive frequency components of 
the spectrum of a pulse of carrier signal of frequency f 0 . 



A 



1/Tp 

-►I 




Fig. 2.8 
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Consider 

i) a sonar operating at a frequency of 3 kHz and R = 15 km 
and ii) a radar operating at 1.5 GHz and R = 150 km. 

The values of pulse repetition frequency are then 1/20 Hz and 500 
kHz respectively. By reference to Fig. 2. 8 it can be seen that in 
sonar the Doppler shift for a modest radial velocity of 1 m/s (4 
Hz) exceeds the line spacing in the spectrum (1/20 Hz) . In the 
case of the radar, the Doppler shift for a mach 2 (600 m/s) 
target is only 60 Hz with a line spacing of 500 kHz. To detect 
the Doppler in the radar requires detection of the shift in the 
spectral line associated wth the carrier. In sonar, it is 
necessary to detect changes in the whole sinX/X spectrum. 

3 ARRAY SIGNAL PROCESSING 
3 . 1 System Model 

A receiving array with a large aperture is directional and 
is able to discriminate against noise sources that fall outside 
the array's directional pattern. An array with an area A operat- 
ing at a wavelength X has a field of view X2 /a steradian. Such an 
array is also beneficial as far as clutter (or reverberation) is 
concerned. Because the array has a restricted field of view, the 
number of clutter (reverberation) scatterers is reduced by com- 
parison to an omnidirectional array whose field of view is 4x 
steradian. Thus, an array is able to improve the signal-to- 
interference level simply by virtue of its size. 

As was stated in Section l.l, an array containing M elements 
suitably spaced is capable of increasing the information rate of 
a system by the factor M. This increase is due to the increase in 
data fed into the system in an M-element array by comparison to a 
single element. If this increased data is processed correctly a 
number of individual beams can be formed each of which is steered 
in a different direction. The M elements sample the spatial field 
and if the spacing is more than half wavelength, the resulting 
beams are said to be independent of each other. This spacing 
requirement is equivalent to the Nyquist sampling theorem. 

Imagine an array in the form of a straight line with the 
elements equally spaced at a distance d apart, see Fig. 3.1. Such 
an array is common in modern sonars where to obtain a sufficient- 
ly large aperture operating at low frequencies a towed array is 
used. In radar also the synthetic aperture that is formed in many 
space-bourne systems forms such a line array. A major difference 
in the processing needs of radar and sonar arrays results again 
from the difference in propagation speeds. 

In sonar it is vital to steer beams because of the otherwise 
slow data rate. For example, consider a sonar with a maximum 
range of 15 km, a sector to interrogate of 360° and resolution of 
10°. It would require 20 s for the sound to travel to maximum 
range and back to the receiver. If the array were to be mechanic- 
ally scanned, the array could not be moved mechanically during 
this time. Thus it would take 36 sequential transmissions of 
sound in each of the 36 bearing resolution cells to interrogate 
the complete 360° sector. The minimum time taken for a single 
look in each bearing cell over all ranges is 720 s or 12 minutes. 

A radar, on the other hand, operating to a range of 150 km 
only requires 10 |is for the electromagnetic energy to travel from 
the transmitter, to maximum range and back to the receiver. Thus, 
for a 1° resolution cell in bearing it is theoretically possible 
to interrogate the complete sector in only 3.6 ms. In practice of 
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Fig. 3.1 

course , the array could not be rotated sufficiently fast and 
simple radars dwell in the same bearing resolution cell for many 
pings of the transmitter. 

Beam steering by electronic means in radar is used to scan 
in elevation while allowing mechanical scanning to provide azim- 
uthal information. The Plessey 3-D AR radar is such a radar. Some 
radars operating at in the h.f. frequency band (where the wave- 
lengths are tens or even hundreds of meters) use electronic 
methods to steer the beams. When such h.f. radars require a high 
angular resolution, it is impossible to rotate the massive arrays 
employed. Sonar must rely upon electronic beam steering in all 
but very short-range situations. 

Suppose that a plain wave approaches the array shown in Fig. 
3.1 from a bearing ©.There will be difference between adjacent 
elements in path lenth travelled by the plain wave equal to 
d.sin©. This results in 

i) a difference in time delay equal to (d.sin©)/c 
or ii) a difference in phase of (27r.sin0) /X. 

The processing necessary for a simple beam steering system is to 
place either the time delays or phase shifts in the appropriate 
signal paths to make all the signals in phase before they are 
added to form a beam output. 

Time-delay beam steering methods are wide-band in their 
operation. This is essential in most sonar operations. Even 
though sonar signals only have bandwidths that are typically say 
5 kHz they are base-band signals and the ratio of centre 
frequency to bandwidth is very low. Time delays in sonar systems 
is readily achieved by using random-access memories that are 
addressed so that they behave like serial shift registers. Some 
radars still use coaxial lines to produce the differential time 
delays. In digital processing schemes, the penalty for using time 
delays to steer beams is that the signals need to be sampled at a 
minimum of 0.5M times the the Nyquist rate. Interpolation tech- 
niques are used to avoid this high sampling rate, but these all 
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use much addition processing power. 

Many phased arrays only produce coherence at a single 
frequency. Hence phasing techniques are generally only narrow 
band. Radars operating with a pulsed-carrier signal of a single 
frequency use phasing techniques to steer beams. One advantage of 
a phased array is that in theory the signals from each channel 
need only be sampled at the Nyquist rate. It should be stressed 
here that the specific case where the phase shift 4> = -toT 
produces a pure time delay Tat all frequencies. However, most 
phasing techniques do not have this property. When such general 
phasing techniques need to be used with broad-band signals, the 
signals need to be separated by filtering into a number of narrow 
bands. Each band is applied to a beam-steering unit where the 
inter-element phase shifts are made proportional to frequency. 

The outputs of these individual beam-steering units are then 
combined to produce the required wide-band output. 

3 . 2 Narrow-band Operations 
3.2.1 Physical interpretation . 

Suppose that a plain wave approaches an array of elements 
from the broadside direction. If the elements form a staight line 
and are equally spaced as in Fig. 3.1, the waveforms produced at 
the output of each element are as shown in Fig. 3.2. If these 
waveforms are sampled at the instant indicated, it can be seen 
that the outputs of the sampling operations produce a signal of 
equal amplitude in each channel. 

Fig. 3.3 shows the effect of sampling the signals received 
on the same multi-element array when the single-frequency input 
waveform produces sinusoidal outputs on the array elements that 
are delayed with respect to each other. It can be seen that if 
these samples are placed so as to form a data series, the series 
is a sampled sine wave. The frequency of the sine wave increases 
as the bearing of the signal increases. As already shown in Fig. 

3.2 a broadside target (bearing - 0) produces samples of a 
zero-frequency signal. Thus, one way of estimating the bearing of 
the target is via spectrum analysis of the sampled data series. 

A problem does arise because the simple process illustrated 
in Fig. 3.2 and Fig. 3.3 produces the same frequency component in 
the sampled data series for targets arising in the port and star- 
board sectors. This left-right ambiguity is overcome by sampling 
the input signals in quadrature. The negative-frequency com- 
ponents represent beams formed in one half sector. Those from the 
other half sector are represented by the positive-frequencies. 

3.2.2 Mathematical representation 

The range of inter-element phase shifts that can be inserted 
in the signal paths before addition to form the individual beams 
is limited to 2n rad. If this value of 2n rad is exceeded, the 
phases repeat themselves because of the periodicity of the sine 
and cosine waveforms involved. With M elements sampling the 
spatial field, no more than M independent beams can be expected. 
If these beams are formed by inserting equal increments of inter- 
element phase shifts, the r th beam will be formed by using an 
inter-element phase shift of 2nr/M rad. The phase shifts used to 
form the r th beam produce a linear phase taper across the array. 
It is sensible to label the array elements starting at one end 0, 

1, 2, ,k , (M-l) and the phase shift applied to the 

k th element is 2Jrrk/M. 
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The r th beam is then formed by summing the phase-shifted 
signals. If a target exists in the direction ©given by the 
equation (2md. sin©) A * 2mrk/M, the phase-shifted signals will add 
coherently to produce a maximum response. This phase shifting and 
subsequent addition may be represented by 



M-l 

B r “Zxk.y.27rrk/M (20) 

k=0 



Note that Xk is a complex number. This was discussed in 
Section 3.2.1. Hence, the spatial samples must be formed by 
quadrature sampling. A more elegant way of representing the phase 
shifts is to use multiplications by cosine and sine in quadrature 
channels. This can be condensed into a single complex exponential 
form. Hence: 

(21) 

k=0 

Of course, this is the same expression as Eqn. 14, the 
discrete Fourier transform. Thus, the process of beam steering in 
narrow-band systems is to form discrete Fourier coefficients via 
the discrete Fourier transform. For this reason, the beams are 
often called spatial frequencies. Algorithms, such as the fast 
Fourier transform, devised to perform the discrete Fourier trans- 
form efficiently and quickly, are equally applicable to both 
sampled time series and to spatial samples. 

4 Correlation Processing 
4 . 1 Averaging 

One method of improving the signal-to-interference ratio is 
to average signals. This assumes that the interference is in- 
coherent and that the averaging will only increase the power 
associated with the interference component in a mean-square 
manner. The signal on the other hand is assumed to be coherent 
from sample to sample. The averaging process will then combine 
the signal components together coherently so that the power is 
proportional to the square of the sum of the magnitudes. Ideally, 
this results in an increase in signal-to-interference ratio prop- 
ortional to the number of samples averaged. 

Correlation processing comes in this general area of time 
averaging. Sonar and radar signals vary with time. For example in 
narrow-band systems the signals can be expected to be sinusoidal. 
Time averaging is only advantageous if the signal being averaged 
is constant from sample to sample. To achieve this the signals 
received in sonar and radar systems must first be modified. The 
method employed is to create a model or replica of the expected 
signal. This is then multiplied by the received signal and the 
averaging process then removes all high frequency components 
leaving an average value of the product with frequencies at or 
very near to zero frequency. This averaged product is the 
cross-correlation between the two waveforms at a particular value 
of time delay. 

Consider an active signal where the transmitted signal is 
y(t) . Let the received signal be x(t). This will be delayed due 
to the time of flight, T, to and from the target. The recieved 
signal is assumed to have the same form as that transmitted. If a 
Doppler shift is expected, y(t) must be modified to account for 
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the Doppler shift. Hence, the model for the process can be the 
transmitted signal or a Doppler-shifted version of this signal. 

To account for the time delay in the received signal, y(t) must 
be delayed by T before being multiplied by x(t) and the product 
averaiged. Hence, we can write the cross correlation coefficient 

o o 

C (T) - J x(t) .y(t+T) .dt (22) 

This represents the correlation between x(t) and y(t) . Strictly 
Eqn. 22 should be divided by a normalising factor to represent a 
true average but this factor is often ignored. In discrete form 
with N samples present and a time delay of d units the discrete 
correlation coefficient becomes 

N-l 

c d = E x k .y k+d (23) 

k=0 

Again Eqn. 23 should be divided by N for it to be a true average. 
This factor is usually neglected. 

Radar signals are often processed using a surface acoustic 
waveform correlator. This correlator uses a replica in sampled- 
data form and a continuous version of the received signal. 
Similarly, sonar systems have used correlators made from charge- 
coupled devices. Here both the signal and the replica are sampled 
data signals but they are not coded in a binary fashion. Neither 
is then of the form normally recognised as digital signal 
processors. However, both sonar and radar regularly use equations 
such as Eqn. 23 to carry out cross correlation using true digtal 
signal-processing techniques. Signals with time-bandwith products 
of the order of 5xl0 3 are quite easily processed. 

4.2 Relationship with the Fourier Transform and Beam Steering 

Comparison of Eqn. 11 with Eqn. 22 and Eqn. 14 with Eqn. 23 
reveals that the Fourier and discrete Fourier transforms are in 
fact correlation processes. The delay terms T and d respectively 
in Eqn. 11 and Eqn. 14 are zero. Thus, these Fourier transforms 
are the cross correlations of the signal with 

(i) a cosine waveform of zero delay in the in-phase channel 
and (ii) a sine waveform of zero delay in the quadrature channel. 
The ratio of the signals in the two channels provides for the 
delay. It is of interest to note that the chirp-z transform [4] 
is an algorithm based upon a sampled data correlator to evaluate 
the discrete Fourier transform. 

Using the same arguements, beam steering is also a correl- 
ation process. Eqn. 21 can be compared with Eqn. 23. It is seen 
that the model assumed for the beam-steering system is that the 
spatial samples will consist of a number of equally-spaced, 
discrete spatial frequencies. The system model is then cross 
correlated with the spatial samples to produce the beams, B r . 

4.3 Correlation Processing in the Frequency Domain 

It is often advantageous to evaluate correlation coeff- 
icients by working with the spectra of the signals. In 
particular, by using efficient algorithms such as the fast 
Fourier transform, the correlation process may be evaluated very 
quickly and with a much reduced computational load. Consider the 
Fourier transform of Eqn. 22 

*o oo 

C(to) = J c(T).exp(-jioT) .dT = J J x(t) .y (t+T) .dt.exp(-jwT) .dT 
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Write (t+T) = s so that ds = dT and separate the variables so that 

oe a0 

C (CO) = J x (t). exp(-jwt) .dt . J y (s) .exp(+jws) .ds 

The first integral is the Fourier transform X(co) of x(t) . The 
second integral is very similar except that the complex expon- 
ential multiplier has a positive sign in the exponent. This 
simply changes the sign of the imaginary part of the Fourier 
transform. Hence the second integral is the complex conjugate of 

the Fourier transform. This is written as Y*(co). Thus, we have 
the simple relationship 

C (co) = X(w) . Y* (w) (24) 

This equation is known as the Weiner-Khinchine algorithm and 
it is used regularly in radar and sonar signal processing to re- 
duce the computational load. The advantage of using the algorithm 
is often missed. In the discrete form, with both data sequences 
of length N, N separate correlation coefficients are evaluated. 
Care must be exercised when using the discrete form of Eqn. (24) , 

C r = X r .Y*, (25) 

since the discrete Fourier transform assumes the data is period- 
ic. Most signals in radar and sonar are aperiodic. To overcome 
this lack of periodicty the data in one of the data sequences can 
have zeros added (zero padding). If this is not done only the 
coefficient with zero time lag is meaningful. When increasing 
numbers of zeros are added to one of the sequences an increasing 
number of coefficients are meaningful. It is usual to take two 
data sequences with N data points each. One consists of 0.5N non- 
zero data samples and 0.5N zeros. The second sequence contains N 
non-zero data points. The application of Eqn. 25 then produces 
0.5N valid coefficients. The remaining 0.5 invalid coefficients 
are discarded. The time windows are made to overlap and in this 
way all the required coefficients are evaluated. 

5 cpndwfl ing Bwmrfcg 

This paper has looked at the differences and similarities 
between the reqirements and execution of digital signal 
processing in radar and sonar systems. The major differences are 
caused by the 2xl0 5 ratio in the velocities of propagation of 
electromagnetics and sound. This results in differences in the 
bandwidths of the signals employed, the durations of the time 
windows that can be processed, differences in the ways that 
Doppler shifts are processed and the urgent need in sonar for the 
use of beam steering using electronic methods. 

In spite of the vast difference in the velocities of 
propagation, surprisingly the wavelengths involved in radar and 
sonar overlap. At h.f. radars also need to employ electronic 
methods to steer beams. At higher frequencies it is becoming 
fashionable to process the signals received by the radar array to 
avoid the necessity of mechanical rotation. The time-bandwidth 
products involved in radar and sonar are also very similar. 

The paper has also reviewed the basic methods of processing 
in the time, frequency and spatial domains. Spectrum analysis is 
often the key to the algoritms used. Where it is relevant a 
physical explaination of the methods described have been given 
together with a brief interpretation of the mathematics involved. 



1/20 




The paper has been purposely restricted to a description of 
techniques based on Fourier methods. More modern methods of 
analysis [5] using alternative models follow in subsequent 
papers. Similarly, details of the methods employed will also be 
covered in other papers. 
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1. INTRODUCTION 

The Digital Signal Processing algorithms employed in radar and sonar applications 
invariably involve the manipulation of time series samples taken from a multiplicity of 
sensor elements. Matrix algebra provides a powerful symbology to represent such 
samples and to analyze their manipulation in a concise manner. Once candidate 
algorithms have been expressed in this algebra there are a wealth of techniques and 
theorems available to assist with the prediction of performance, and thereby optimise the 
design of the finally implemented signal processing structure. 

In what follows we review the use of matrix algebra to represent sampled signals 
and to analyze Digital Signal processing algorithms. We also present some of the key 
theorems of matrix algebra that find applications in such analysis. 

2. THE VECTOR REPRESENTATION OF SAMPLED SIGNALS 

Consider an arbitrary band-limited waveform, x(t), originating from a single sensor 
element. If v ^ denotes the maximum frequency component in x(t), it follows from the 
Nyquist sampling theorem that the waveform can be re-constructed for all values of t 
from samples taken at times t k = k*t, with integral k and At < 1/2^. If the bandwidth 
of the waveform, B, is less than one may seek to reduce this sampling interval by a 
suitable down-conversion process. In the limit one can in principle down convert such 
that the waveform spans the baseband range o^u^B. In this case the Nyquist sampling 
interval is reduced to At = 1/2B. & 

A series of p such samples with m-p+l<k£m may be represented by a p 
component column vector, x given by 

* " ( X m-p+ V X m-p+2 X m) T 

where x k is the sample obtained at time kAt. The notation T denotes the transpose of a 
matrix. For an arbitrary matrix this operation involves the interchange of rows and 
columns. It is used here because it is more compact to represent vectors in row form. 

An alternative to using real samples is to represent the waveform in complex 
notation. The waveform x(t) can be viewed as being composed of sinusoidal components 
and expressed in the form 

*STC Technology Limited, London Road, Harlow, Essex, CM17 9NA 
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x(f)- J"i4(w)cos[ut+4>(w)]Ju. 



We may then define the complex representation of the waveform 
by 



x(t) - J*i4(w).exp{ j [wt+4>(w)] }<fu. 



We then have that 



x(t) - R e { x(t) }. 



We also have that 



I m { x(t) )- j m A(w) sin [wr+<J>(w)]du. 

- f~A(w) cos [wt+4>(w)-7i/ 2 ]du. 



Samples of the imaginary component of X(t) can therefore in principle be 
obtained from x(t) by passing the latter waveform through a filter with a flat unit 
amplitude response and a frequency-independent phase shift of -rr/ 2 . Alternatively 
samples of the imaginary component can be directly computed from samples of x(t) by 
means of the Discrete Hilbert Transform which performs this filtering operation digitally. 

We may therefore alternatively represent the series of p samples with 
m-p+l<k<m by a complex column vector given by 

i J > *-,.2 ••• *J T - 



A third representation that is commonly employed involves the use of complex 
baseband notation. The complex baseband representation of the waveform x(t) may be 
defined by 
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*(0 - *(<)• exp i-j w/}. 



Here w Q is the angular frequency at the centre of the band of interest. 
fi(t) may be expressed in the form 

x(t)- J* m A(w) exp {jwtidv, 



where A(w) is a complex function of w spanning both positive and negative frequencies. 

The real and imaginary components of tf(t) have spectra which span the band 
o<v<B/2. The Nyquist sampling interval appropriate to each is therefore At = 1/B. This 
is twice the interval required for real baseband samples, but each complex baseband 
sample has both a real and imaginary component. 

The real and imaginary components of £(t) can be obtained by mixing x(t) with 
orthogonal local oscillator signals of the form cos[w G t] and -sin[w 0 t] and subjecting the 
resulting waveforms to low pass filters which span the frequency range o<u<B/2. This 
mixing and filtering may be undertaken in the analogue domain, in which case 
appropriate samples need to be taken at the outputs from the low pass filters, or in the 
digital domain in which case the input to the process will be the sequence x k of real 
sample values of x(t). 

We have thus far developed three possible vector representations of the last p 
samples obtained from a single element sensor. If we have n such sensor elements we 
may similarly represent the set of samples obtained at time t = mAt across the array by 
an n component column vector of the form 

* - (*!,.« * 2 * - X *J T ' 



Here x^ is the sampled output from the ith sensor at time t=mAt. Such x^ can 
either be real samples, or samples in one of the two forms of complex representation we 
have considered. It is now a simple extrapolation to represent the last p samples from 
all the n sensors by an np component column vector. 

3. DIGITAL FILTERING IN VECTOR NOTATION 

Now that we have a column vector representation of signals sampled both in time 
and, via an array of sensor elements, in space, we may consider the development of 
sampled frequency and spatial filters in vector notation. 
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Consider the sampled output from the ith sensor. We may apply a sampled finite 
impulse response filter to the output from this sensor by forming a weighted combination 
of the last p samples. We then have for the output, y, appropriate to time t=mAt 

y- E w m-i w 

j-o 



Here the Wj are the components of the impulse response of the filter. We may represent 
this impulse response by a p component column vector, w, given by 

“ - K-,.p w -) r - 



We then have in vector notation 



y - 34! r £. 



where 



£ 



(*i 



Um-p+V ijm-p+l* 






If we are working in one of the two complex notations, the components of w will 
be complex, representing amplitude scaling factors and phase shifts applied to the 
complex representations of the sensor samples. Conversely, if we are working in real 
notation the components of w will be real. The vector notation for the filter output is 
equally valid in all three cases. From now on we will develop the matrix algebra for the 
more general complex case. The processing of real samples can then be treated as an 
appropriate special case. 

The finite impulse response filter we have applied to the output from the ith 
sensor implements a filtering operation in the frequency domain. If we replace x by the 
vector of samples obtained across the array at time t=mAt we may similarly implement 
a spatial filter. Such a filter corresponds to a far field reception pattern for the sensor 
array which varies with bearing. In general, with a column vector, x, containing both 
space and time samples we may implement a combination of frequency and spatial 
filtering. In all cases we may employ the simple inner vector product 



y - & T & 



to represent the filtering operation. 
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4. FREQUENCY AND SPATIAL FILTERING 

When the components of the column vector x represent the last p samples taken 
from a sensor element, the transformation y = w T x corresponds to that of implementing 
a discrete time finite impulse response filter in the frequency domain. If the samples are 
in complex baseband notation, a particularly important form for w is one of the set y^, 
with o<q<p-l, given by 

yv - — (1, exp{-;— q\ exp i-j— (p-l)#}) 7 . 

q y/p P P P 



Such an impulse response seeks to combine coherently the time series samples of 
any component of the sensor waveform that has a frequency given by 

q k 

l) - -2_+ — 
ptX At 

where At is the sampling interval, and k is an appropriate integer selected such that 

1 1 

+ — . 

2 a t 2 a t 



For a more general component of the form S c exp{jwt}, the output from the filter 
at time t= kAt follows as 






io 

€ 



sin 

sin [n(qlp - da/)] 



The complete set of Wq thus corresponds to a set of filters with centre frequencies 
separated by au - 1/pAt which span the band 



1 

2aT 



+■ 



2aT 



If we define 



E - te e , is;,,... 
and 

x - fy,. V- V.) r * 
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we may write 



y - £ i. 



The column vector y then has components equal to the outputs from each of the 
filters in the set. This equation expresses, in matrix notation, the operation of taking the 
Discrete Fourier Transform of the samples that are the components of x. It is 
straightforward to show that F is a full rank transformation. It therefore has an inverse, 
F 1 , defined by F 1 F = I where I denotes the identity matrix. We may therefore write 

* - ET 1 y. 



It follows from Fourier Transform theory that if we set 
f - iv, ... * 

then 

*> 'jf ('• “ p Hr '}• Ff 2 '}- FfH)' 

That is 



F 1 - F. 



If the components of x correspond to a single time snapshot of samples from an 
array of sensors, the transformation y = w T x corresponds to that of implementing a 
spatial beamforming operation. The spatial domain equivalent to a single frequency 
component is a plane wave incident on the array from a specific bearing. Let such a 
wave have an angular frequency w. The output from the ith sensor may then be 
expressed in the form a x exp {jw[t-rj}. If the receiver chain associated with each sensor 
produces noise of equal power, but uncorrelated from sensor to sensor, it follows from 
matched filter theory that the optimum weight vector to combine these outputs is given 
by 



yy - “H («i exp { j w Tjl, a 2 exp {j w t 2 }, .... a n exp {j w t b }) r , 

where 17 is an arbitrary scaling factor. For the case of a linear array of equally spaced 
isotropic elements the Discrete Fourier Transform forms a useful set of receiving beams 
spanning all space. 
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5. THE REPRESENTATION OF POWER THROUGH COVARIANCE 
MATRICES 

The instantaneous value of the magnitude squared of the voltage at the output of 
an arbitrary filter follows as 

y* y - (i v T (tt r *). 

Here the notation * denotes the operation of complex conjugate. 

Since w T x is, in general, a complex scalar, it follows that 

w T x - (& T x) T . 

Now, for arbitrary matrices A and B it may readily be shown that 

fi r A r . 

Thus 

yy T x - (& T x) T - x T )&>. 



We thus have that 

y*y - (}v T £ r iv - yz H x* & T iv. 

Here the notation H denotes the operation of taking the complex conjugate of the 
transpose of a matrix. 

The expectation of the magnitude squared of the output voltage therefore follows 
as 

P - E { y* y } - w. H E { x* x T } m 



Here E{...} denotes the expectation operator, and it has been assumed that w is time 
invariant. We may define the covariance matrix, R, characterising the components of x 
by 



R - E { x* i r }. 



We then have that 

E H - E { (x* xY ) - E { (x x B ) T } - E { i* x T ) - & 
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Such a matrix is said to be Hermitian, or in the case of real variables orthogonal. 
In terms of the co-variance matrix we thus have 

P - w. B R w. 



P is termed the filter output power. In the case of real variables this is an exact 
definition if the sampled voltages are referenced to a transmission line of unit 
characteristic impedance. In the case of complex notation, P is strictly twice the power, 
but this factor of two is conventionally ignored. 

The diagonal components of the co-variance matrix R represent the powers in the 
individual components of x. The off-diagonal elements represent correlations between 
distinct components of x. 

6. CONSTRAINED MINIMISATION IN MATRIX ALGEBRA 

A powerful tool much used in the solution of constrained minimisation problems 
with real functions of several variables is that of the method of Lagrange Multipliers. If, 
for example, we wish to compute the values of x and y that minimise 0 = f(x,y) subject 
to the real constraint y=y(x), rather than eliminating one of the variables in f(x,y) with 
the aid of the constraint equation, we can instead consider the function 

® - / too + * t y -y <*)L 



where A is a Lagrange multiplier. A necessary condition for the constrained minimisation 
of 0 is that 



do 

dx 




These two equations, together with that defining the constraint, can be used to 
solve for the three unknowns x,y and A. The procedure can be extended to deal with 
multiple constraints by adding additional Lagrange Multipliers. Each additional 
constraint introduces an additional unknown Lagrange Multiplier and an additional 
constraint equation such that the number of equations remains equal to the number of 
unknowns. 

Quantities to be minimised are necessarily real, but constraints are in general 
complex. Suppose for example we wish to find the weight vector, w, that minimises the 
output power y^Rw, from a filter, subject to the constraint that a certain weighted 
sum of elements of the weight vector is equal to a given, in general complex, number. 
Such a constraint may be written in the form 
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E c . w rP- 

i 



If we define 

the constraint equation can be expressed in matrix form as 

£ r lv-P. 



This complex constraint can be viewed as two real constraints, given by 

R e 



and 

I m k T Yd-I m < P>. 

We are therefore led to consider the function 

0 - w B R ie+Xj R e k T ii?-p> + X 2 l m k T -p), 

where X 1 and X 2 are two real undetermined Lagrange Multipliers. This augmented 
equation can be re-expressed in the form 

0 - UJ* £ W. + — \sJ 14^-p + £ H iH*-P*] 

M>-P -c B 



w" R w. + t 



*1 



+[■ 



-j y] [ q t k-P] 
V y] k B 1K*-P*] 
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Hence 

0 - *!* R M!+X[ C T M!-P] + A.*[c" Mf-P*]- 



Here 



* 




-j 



hi 

2 



is a single complex Lagrange Multiplier. 

In order to proceed to obtain the solution to the constrained power minimisation 
problem we must now evaluate the complex vector w that minimises 0. 

We may write 

0 - 



where 






and 



- y±>j - j & Q 



Thus 



& - & te, *j W.J -j YV.J 



Now 



* 0-E 



Jh 

dW r 



6 w. 



Jl 

dw 



6 w 



(?< 
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where w u is the ith component of w I? and the ith component of Wq. 

We may now define a complex vector gradient, v 0, by 

- * f dt] . dn dn . dt] V 

* 0 ^ ^qi ty n ^ ^02 ] * 



which is a generalisation of the three dimensional real gradient operator denoted by grad 
or V. For a small change 5w = 5wj +j 5 wq in the weight vector w the corresponding 
changes in 0 follows as 

60 - R t {(iz o) r ) 6 + I M (Cy e) r ) 6 w. Q 

- \ [Ci0) T + Cia)'] & VLj - \ j KS 0) T - (5 a)®] 6 T6L q 

- -i {(s e) T [8 M! ; - j 8 + (20)® [8a! / * j b 

- R e ((50 ) T 6 - R e 0} 



It therefore follows that, if v 0 = <^50 = 0 for all 0 w, the condition required 
to find the desired minim um of 0. 

Since 0 is expressed as a function of w and w\ rather than w, and Wq, it is 
convenient to re-express v 0 in terms of derivatives with respect to the components of 
w and w\ 

Now 

dn da ^ da 9 wi 
dw B dw ( dw B + d>v ; dw B 
da da 
dw t dw’' 

We also have that 



an da dw, t da dw- 
dw 



. 00 00 
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Hence 



0ri .dii - 0e> 
C20) ( - -r- 1 - + J T 1 - " 2 

Qi dw’ 



Thus 



^ 0-2 



^ 00 00 



\T 



0 ^; 0 *; j 



Returning to the expression for 0, we have that 

0 - yv H R iv+A (c r it -P) + r(c* iv*-p*). 



Hence 



vo - 2 [fi iv + r £* ]. 



The required condition for the minimisation is therefore 

ve> - 2 [fi tt+r £ * ] - a 



Hence we require 



R YL - -r £*. 



Pre-multiplying both sides by the inverse of R, we thus have that 

yv - -V RT 1 £*. 



But 



c T w - -X* £ r R- ! c* - p. 



Hence 



A* 



C r S' 1 c* 
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Thus 



w. - 0 



ir 1 g* 
C T BT l c • 



7. EIGENVALUES AND EIGENVECTORS 

Let us consider a further constrained minimisation problem, namely that of finding 
the weight vector w that minimises the output power subject to the sum of the squares 
of the moduli of the components of w being constrained to unity. This constraint can 
be expressed as 



W. H w. -1 



and is necessarily real. We therefore consider the real function 

0 - iv* R -A, (& H iv -1), 

where is a real Lagrange multiplier. We thus require 

- 2 [ Rw - - Q, 



or 



R w. - A,*. 



We are thus seeking a column vector w which when pre-multiplied by the matrix 
R produces a scaled version of w. Such a vector is called an eigenvector of the matrix 
R, and the scaling factor is termed the corresponding eigenvalue. In general there are a 
number of eigenvectors of a matrix. Let the ith be denoted by Wj and its corresponding 
eigenvalue by We then have that 

* 



Clearly, if Wj is as eigenvector of R so is any scaled version of Wj. We may 
therefore, without loss of generality, choose to work with normalised eigenvectors 
satisfying Wj H Wj = 1. 

Pre-multiplying both sides of the defining equation for Wj by w” we find that 
ivf R jv. - A., jvf iv.. 
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Taking the complex conjugate of the transpose of both sides we then have that 



xl b R H yv ( - X* yl b w. c 



When, as in the case of a covariance matrix, the matrix R is Hermitian we have that 



R H - & 



Hence 



yvf R XL. - X] w* yv*. 

We therefore conclude that X* = The eigenvalues of an Hermitian matrix are 
therefore always real. 

Let be a second eigenvector of R, with eigenvalue Xj. Then 

R " X J 



Thus 



XL B R XL j - Xj yl* Wy 



But 

wf R Wj - (R xl) b xl } - ( X { xl^ b ^ - X] xl b xl. - X, xl b xl. 

Hence 

x i < < */ 



Thus, if \ * Xj, we must have that 

xl b xl. - 0 . 

Thus eigenvectors of an Hermitian matrix corresponding to distinct eigenvalues 
satisfy Wj H ^ = O. By analogy with the expression for the dot product of two real 
vectors, eigenvectors satisfying such a condition are said to be orthogonal. 
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Even if the eigenvectors have the same eigenvalue, in which case the eigenvalue 
is said to be degenerate, we may still incorporate them into a set of orthogonal 
eigenvectors. As an example let Wj and Wj both have the same eigenvalue A. Then a Wj 
+ is also an eigenvector with eigenvalue A. We may take as our two eigenvectors 



yv i and a + P 



we may require that 



yvf (aiv j + P 1 1 ^) - O (< orthogonality ) 



and 



(« + p ) B (a + P yn ) - 1 (normalization). 



If we solve these two expressions for a and p we will have constructed two 
orthogonal and normalized eigenvectors. We may generalize this approach to greater 
than two-fold degeneracy to prove that, even in the case of degeneracy, the eigenvectors 
of an Hermitian matrix can be constructed to be normalized and orthogonal to each 
other. They are then said to constitute an orthonormal set. 

We may solve for the eigenvalues of a matrix R by means of the following 
procedure. Since 






we have that 

u- 



where I is the identity matrix with elements 5^. 

If [R - AI]" 1 exists, we have that 

^ - b - * z r 1 a - a 



In this case the only solution for w* is the trivial solution Wj=Q. We may therefore 
demand that [R - A I]' 1 does not exist to gain information about the non-trivial solutions. 
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But 



LB - X zr 1 



adj m - x n 

det m - \ H 



Here adj[ — ] denotes the adjoint of a matrix, and det [...] the determinant of a matrix. 
Since the adjoint of a matrix may always be evaluated, we conclude that the condition 
for the existence of non-trivial solutions reduces to 

det LB - X fl - O. 

The expansion of this determinant when R is an (n x n) matrix will yield an nth order 
polynomial. There are therefore in general n eigenvalues, and hence n eigenvectors, of 
such a matrix. As we have seen above, we may construct these eigenvectors to form an 
orthonormal set. 

The n orthonormal eigenvectors of the matrix R may then be written as the rows 
of a matrix B. We will then have that 

RR T - R T H 



where D is a diagonal matrix containing the n eigenvalues X r 
Since the rows of B form an orthonormal set, it follows that 

R' R T - L 



and hence that 



RR H -l 



A matrix B satisfying such a condition is said to be a unitary matrix. The inverse of a 
unitary matrix is therefore equal to the complex conjugate of its tranpose. 

It therefore follows that 



R' RR r -R' R T n-& 



A transformation of a matrix R by means of a unitary matrix B as in this last 
equation is called a unitary transformation. We have thus shown that any Hermitian 
matrix can be transformed into a purely diagonal form by means of a unitary 
transformation. Such diagonalization of an Hermitian matrix is a powerful tool which can 
be used to uncouple matrix equations as we shall now see. 
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8. THE SOLUTION OF TIME DEPENDENT COUPLED EQUATIONS 



In Section 6 we considered the determination of the weight vector corresponding 
to the minimum output power subject to the constraint c T w = /3. One method of 
evaluating such a vector is to implement a recursive procedure in which the gradient of 
the output power with respect to the currently applied weight vector is computed, and 
a small step applied in the direction of steepest descent until the minimum condition is 
reached. Clearly such steps must not be allowed to violate the constraint condition.The 
gradient can, for example, be evaluated by monitoring the effects of small perturbations 
in the applied weight vector components. 

Such a recursive update procedure may be expressed in the form 

- & k - \i E 1 fef £ n?p. 



Here 



E - 



l- 



£ £ 
£ T £ 



Now 



C T £ A WrC. T 



Z 




A n? 



Hence 

£ T E An? - - c 1 ) An? - Q. 



The matrix P ensures that the desired constraint is maintained by such a procedure 
provided that Wo satisfies c T ^ = 0. P is known as a projection matrix. 

We therefore have that 

n^ +1 - w k - 2p £ R & k . 



The recursive procedure will reach a stead-state condition in which 

* ** “ *- 
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We therefore have that 



E R m. - Q- 

m 



We may therefore write that 

^ +1 - W m - YL k " YL m - 2 p £ R [ YL k ~ JKj. 



If we define 



A yv t - w. k - 



the equation can be expressed in the form. 

A m^ +1 - A w. k - 2 p E R A jv^. 



Since 



C r - C r - P 



It follows that 



q t \& k - uj - C r Aj^ - 0 



Hence 



EAw k - Au^. 



We therefore have that 



A - Am^ - 2 p £ £ £ A iK t . 



This represents a set of coupled equations. We may seek to uncouple these 
equations by noting that P R P is an Hermitian matrix. If B is a matrix with rows equal 
to n orthonormal eigenvectors of P R P we then have that 

&' (ERE) i T - & 
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where D is a diagonal matrix with elements equal to the eigenvalues of P R P. 

Premultiplying both sides of the coupled equation set by 8\ we find that 

fi* A A m k - 2 p g (E R O A yv k 

-f Att r 2|i|*(££Sf (I* A 

- R’ AKj 

since 

& T - CfiV. 



We now have a set of uncoupled equations which we can solve for 8* A w k , and 
hence for A w k . The time constants of convergence of this equation set are directly 
related to the eigenvalues of the matrix PRP. 

9. LEAST SQUARES PROCESSING 

Suppose we wish to find the weight vector which minimizes the output power from 
an arbitrary spatial/frequency filter, subject to the constraint of one of the components 
of the weight vector being held fixed at unity. Although we may address this problem 
using the techniques described in Section 6, it is instructive to re-couch it as follows. 

Let successive snapshots across the filter inputs to which the variable weight 
components are to be applied be denoted by column vectors £. Let the corresponding 
scalar snapshots associated with the input to which the fixed weight is to be applied be 
denoted by y^ Let N such snapshots be available, with l<j<N. We may then define 

and 

i - ( y v y 2 y K ) T 



Let w denote the vector of variable weight components we are seeking. If we 
consider 

e - X w. + y. 
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the best estimate of w after the N snapshots is that which minimizes e H e. The problem 
has therefore been reduced to that of classical least squares processing. 

Now 



£ B £-(Xx + x) B (Xll! + X) 

- * B X B X k + X H y + X* X * + y* y. 



Hence 



V (£ B £) - 2 [X B X *! + X B JO- 



The required value of w therefore follows as 

* - -CX" XT' (X* *)■ 



A problem with using this expression to evaluate w is that the equation set may 
be ill conditioned. That is, small perturbations applied to the input data may produce 
large fluctuations in the solution. Such perturbations arise from the finite step size 
associated with the analogue to digital conversion process. The problem is made worse 
by^the particular equation set we are considering because it involves the computation of 
X X> which means that the input data samples are effectively subjected to a squaring 
process. 

Instead of employing such a squaring process, we can instead consider multiplying 
both sides of the equation defining e by a series of unitary matrix transformations. Let 
the first of this series be denoted by |Q. Since Q H 0 = I for a unitary matrix, we have 
that 



The length of the vector 0 e is therefore the same as that of e. By means of an 
appropriate choice of the series of rotations, the compound effect of which we denote 
by the matrix Q, we may reduce the equation set to the form 

Qi - QXxi + Qx 
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Here R now denotes a matrix with zeros below the main diagonal. It is clear that the 
condition corresponding to the minimum length of Q e, and hence of e, is 

R iv - -b v 



The vector w can then be found by the process of back substitution. This method 
of solving the least squares problem is known as "QR processing". 

A particularly important form for the constituent rotations 0 is when each is 
constrained to operate on only a pair of rows of the equation set at a time, leaving all 
the other rows unaltered. We may represent each such a transformation by a [2x2] 
matrix, which operates on its associated pair of rows, and is given by 




5 

C 



i 



We require that £ H £ = L and hence that 



s*s + c*c - 1 



and 



c* - c. 



These two conditions will be satisfied if we select 



c 



vfl x P + I y P 



and 



s 



y 

X 



c. 



With these forms for c and s we also have that 
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A series of such Givens' rotations may therefore be used to annihilate successive 
components of the data matrix X, until it is reduced to the desired upper triangular form. 

By avoiding the explicit computation of X H X the data domain QR processing 
technique offers a numerically superior algorithm for determining the required weight 
vector. 



10. CONCLUDING COMMENTS 

In the preceding sections we have given an introduction to the use of matrix 
techniques in representing and analyzing time series samples taken from the elements of 
a sensor array. Matrix representation provides a concise notation with which to represent 
the samples taken with such systems, and the wealth of available theorems relating to 
manipulations in this algebra means that it also represents an extremely powerful 
analytical tool. 
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